/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		segment_shape_function.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-22 22:54
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:38
*
*   Description:
*
=============================================================================*/
#ifndef SEGMENT_SHAPE_FUNCTION_H
#define SEGMENT_SHAPE_FUNCTION_H

#include "../../types.hpp"

namespace geoxfem {
namespace geometry {

template <size_t dim, size_t degree>
struct segment_shape_function {
  static constexpr _real r_5 = .5;
  template <typename Shape>
  static VectorXR getStdNByLocal(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
  template <typename Shape>
  static MatrixXR getStddNdlxByLocal(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
  template <typename Shape>
  static VectorXR getStdNByGlobal(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
  template <typename Shape>
  static VectorXR global2local(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
};
/*segment| dimension: 2D, order: 1{{{*/
template <>
template <typename Shape>
VectorXR segment_shape_function<2, 1>::getStdNByLocal(const VectorXR& lcoord,
                                                      const Shape&) {
  Vector2R    N;
  const _real xi = lcoord(0);
  N << r_5 + xi / r2, r_5 - xi / r2;
  return N;
}
template <>
template <typename Shape>
MatrixXR segment_shape_function<2, 1>::getStddNdlxByLocal(const VectorXR&,
                                                          const Shape&) {
  MatrixXR dNdlx(2, 2);
  dNdlx << r_5, r_5,  //
    r0, r0;
  return dNdlx;
}
template <>
template <typename Shape>
VectorXR segment_shape_function<2, 1>::getStdNByGlobal(const VectorXR& gcoord,
                                                       const Shape&    shp) {
  Matrix2R m;
  Vector2R v;
  _real    x1 = shp.front().getCoord().sum();
  _real    x2 = shp.back().getCoord().sum();
  m << 1, 1,  //
    x1, x2;
  v << 1, gcoord.sum();
  return m.inverse() * v;
}
template <>
template <typename Shape>
VectorXR segment_shape_function<2, 1>::global2local(const VectorXR& gcoord,
                                                    const Shape&    shp) {
  Vector2R lcoord;
  Matrix2R lcoords;
  lcoords << -1, 1,  //
    0, 0;
  lcoord = lcoords * getStdNByGlobal(gcoord, shp);
  return lcoord;
}

/*}}}*/
/*segment| dimension: 2D, order: 2{{{*/
template <>
template <typename Shape>
VectorXR segment_shape_function<2, 2>::getStdNByLocal(const VectorXR& lcoord,
                                                      const Shape&) {
  _real xi = lcoord(0);
  _real t1 = -1. - xi;
  _real t2 = -xi;
  _real t3 = 1. - xi;

  Vector3R N;
  N << t2 * t3 / 2.,  //
    -t1 * t3,         //
    t1 * t2 / 2.;

  return N;
}
template <>
template <typename Shape>
MatrixXR
segment_shape_function<2, 2>::getStddNdlxByLocal(const VectorXR& lcoord,
                                                 const Shape&) {
  _real xi = lcoord(0);
  _real t1 = -1. - xi;
  _real t2 = -xi;
  _real t3 = 1. - xi;

  MatrixXR dNdlx(2, 3);

  dNdlx << -(t3 + t2) / 2., (t3 + t1), -(t2 + t1) / 2.,  //
    0, 0, 0;

  return dNdlx;
}
template <>
template <typename Shape>
VectorXR segment_shape_function<2, 2>::getStdNByGlobal(const VectorXR& gcoord,
                                                       const Shape&    shp) {
  Matrix3R m;
  Vector3R v;
  _real    x1 = shp.getPoint(0).getCoord().sum();
  _real    x2 = shp.getPoint(1).getCoord().sum();
  _real    x3 = shp.getPoint(2).getCoord().sum();
  m << 1, 1, 1,  //
    x1, x2, x3,  //
    x1 * x1, x2 * x2, x3 * x3;
  v << 1, gcoord(1), gcoord(1) * gcoord(1);
  return m.inverse() * v;
}
template <>
template <typename Shape>
VectorXR segment_shape_function<2, 2>::global2local(const VectorXR& gcoord,
                                                    const Shape&    shp) {
  Vector2R lcoord;
  Matrix2R lcoords;
  lcoords << -1, 0, 1,  //
    0, 0, 0;
  lcoord = lcoords * getStdNByGlobal(gcoord, shp);
  return lcoord;
}

/*}}}*/
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
}  // namespace geometry

}  // namespace geoxfem

#endif /* SEGMENT_SHAPE_FUNCTION_H */
